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ABSTRACT 

This study, using publicly available simulations, focuses on the characterisation of the 
non-Gaussianity produced by radio point sources and by infrared (IR) sources in the 
frequency range of the cosmic microwave background from 30 to 350 GHz. 

We propose a simple prescription to infer the angular bispectrum from the power 
spectrum of point sources considering independent populations of sources, with or 
without clustering. We test the accuracy of our prediction using publicly available 
all-sky simulations of radio and IR sources and find very good agreement. 

We further characterise the configuration dependence and the frequency behaviour 
of the IR and radio bispectra. We show that the IR angular bispectrum peaks for 
squeezed triangles and that the clustering of IR sources enhances the bispectrum 
values by several orders of magnitude at scales i ^ 100. At 150 GHz the bispectrum 
of IR sources starts to dominate that of radio sources on large angular scales, and it 
dominates over the whole multipole range at 350 GHz. 

Finally, we compute the bias on /nl induced by radio and IR sources. We show 
that the positive bias induced by radio sources is significantly reduced by masking the 
sources. We also show, for the first time, that the form of the IR bispectrum mimics 
a primordial 'local' bispectrum /nl- The IR sources produce a negative bias which 
becomes important for Planck-like resolution and at high frequencies (A/nl ^ —6 
at 277 GHz and A/nl ^ -60-70 at 350 GHz). Most of the signal being due to the 
clustering of faint IR sources, the bias A/^^ is not reduced by masking sources above 
a flux limit and may, in some cases, even be increased due to the reduction of the 
shot-noise term. 
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1 INTRODUCTION 

In the last few decades, the Cosmic Microwave Background 
(CMB) has become a very successful probe of the early and 
late universe. The smallness of the perturbations in the cos- 
mic fluids, and hence in the space-time metric, allows us to 
use linear perturbation theory to compute their evolution ef- 



all these measurements, combined with probes of the geom- 
etry of the Universe , like Baryonic Acoustic Oscillations 
(e.. 



Percival et al. (20101, Blake et al 



Supernovae (e.g. Astier et al. (20061, Hicken et al 



( 2011[)), t ype l a 
(|2009|, 



Guy et ar](|2010[)) and Hubble constant measurements (e.g 



Since COBE (Smoot et al.|1992l the measurement of 



ficiently and accurately through Boltzmann codes ( Seljak & 
|Zaldarriaga] ( [T996l ) [Lewis et al.|(f2000l) [Les gourgucs (2011)) 



Freedman et al. (20011, Riess et al. ( 2009 1 , .Frccdman et al 



'11 



(20091), give a converging view of our Universe and have 



lead to the establishment of a 'standard model' of cosmology 
(e.g. Larson et al. (2011[)) known as ACDM. In this model, 



the CMB power spectrum has been achieved by many ex- 
periments and over a wide range of scales. The most recent 
CMB data are those of the Seven- Year Wilkinson Microwave 
Anisotropy Probe (WMAP) ( [Larson et"ar]|2011[ ), Atacama 



Cosmology Telescope (ACT) (Das et al. 20111 and South 



Pole Telescope (SPT) ( [Keisler et al.|2011[ ). Constraints from 



the universe is flat, dominated by a cold dark matter com- 
ponent and a 'dark energy' component compatible with a 
cosmological constant. The present constraints suggest that 
CMB anisotropics are a realisation of a primordial random 
process that generated the initial perturbations from quan- 
tum fluctuations which were then str e tched t o macroscopi c 
scales by inflation (e.g. 



(19791, 



Guth 



(19811, 



E-mail: Fabien.Lacasa@ias.u-psud.fr 



Starobinskii 
Liddle fc L^ ( [2000[ ), see [Bassett et al.[ ( |2"o66| ) and ^Linde] 
(2008[) for reviews). 
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The microwave sky is however not made of CMB pri- 
mary signal alone, it consists also of secondary anisotropies 
such as those associated with the Integrated Sachs- Wolfe ef- 
fect - arising from the time-variable metric perturbations 
(Sachs & Wolfe 19671 -, those arising from the Sunyaev- 



(20031, Creminelli fc Zaldarriaga (20041) that is sub- 



Zel'dovich (SZ) effect (inverse Compton scattering) in the 
direction of galaxy clusters (Sunyaev & Zeldovich 19721, 



those due to the Doppler effect from moving structures (e.g. 
kinetic SZ effect from clusters and inhomogeneous reionisa- 
tion), see Aghanim et al. (20081 for a review on secondary 



anisotropies. In addition, there are other astrophysical com- 
ponents in the microwave domain such as the dust, syn- 



chrotron and free-free emissions from our Galaxy (Planck- 
Collaboration|[2011e I, the emission from radio sources that 



dominate at lower frequencies but contribute significantly at 
microwave frequencies ( de Zotti et al.|[2005 l, and the emis- 
sion from dusty star-forming galaxies emitting mainly in the 
Infra- Red (IR) domain ( [Low fc Tucker|1968[ ). 

In this study, we will focus only on the characterisa- 
tion of the extra-galactic point sources, namely the radio 
sources and the IR dusty galaxies. They contribute notably 
to the power spectrum at CMB frequencies and start domi- 
nating over the CMB at about £ ~ 2000 when the CMB sig- 
nal plummets. Active Galactic Nuclei (AGN) are observed 
as radio sources via their synchrotron emission. They have 
been widely studied in the CMB context especially at low 



frequencies !/ ^ 90 GHz (e.g. de Zotti et al. (20051, Boughn 



fc Partridge (20081, Sajina et al. (2011 1). They affect mostly 



the lower end of frequencies observed by CMB experiments. 
Their largest impact was thought to be in the frequency 
bands from 30 to 90 GHz but the recent Planck results 
( Planck-Collaboration|2011b|a I show that radio sources are 
detected at frequencies as high as 217 GHz. At the CMB 



frequencies, radio sources do not cluster (Toffolatti et al 



( 1998 1 Gonzalez-Nuevo et al. ( 2005 l) and thus exhibit a flat 
power-spectrum (see Appendix |A[). 

Star-forming galaxies are observed as IR sources via the 
thermal emission from dust heated by the ultra-violet emis- 
sion of young stars. These IR sources contribute to the signal 
observed by CMB experiments at frequencies higher than 
150 GHz, thus they are particularly relevant for the most 
recent CMB experiments, e.g. SPT, ACT and Planck High 
Frequency Instrument observations. The cumulated emis- 
sion from the IR sources, the Cosmic Infrared Background 
(GIB), was first discovered by _Puget et al. (1996J, and its 



anisotropies were first characterised by [Lagache fc Puget| 



(20001 and Matsuhara et al. (20001. Many other observa- 



tions were possible in the last decade in the IR and submm 



et al. (20101, Amblard et al. (20111, Planck-Collaboration 



domain (Lagache et al. (20071, Viero et al. (20091, Hall 



(2011dl). In particular, the latest constraints of the CIB 
from Planck showed that its power spectrum, at frequen- 
cies 217-353-545-857 GHz, is well modeled by a power law 
C7f «(i.) = A(u)x{j^^)" with e.g. A = (104±4)xl02jyVsr 
and n — —1.08 ±0.06 at 545 GHz. This behaviour contrasts 
with the flat spectrum of radio sources and is due to the 
clustering of the IR galaxies and their host dark-matter ha- 
los. 

These point sources are super-imposed on the pri- 
mordial fluctuations. The simplest models of inflation 
(single-field and slow-roll) predict a small primordial non- 
Gaussianity (NG) (Acquaviva et al. (20031, Maldacena 



dominant to the NG induced by the non-linear evolution of 
the perturbations, a contribution that is necessarily always 
present. More complex inflationary models, e.g. multi-field 
scenarios, may predict larger NG ( Byrnes fc Choi||2010 l, to 
the point of being detectable. A simple and yet powerful 
probe of non-Gaussianity is the three-point function in har- 
monic space, the angular bispectrum (see sectionplfor more 
details) , which is defined as a function of a multipole triplet 
(^i,^2,i'3). The angular bispectrum vanishes for a Gaussian 
field like all odd-order moments. Besides the bispectrum, 
connected even-order moments may also be used to probe 
non-Gaussianity, and the 4-point function or trispectrum 
has indeed been a focus of interest ( Kunz et al.|2001 l, espe- 
cially for lensing studies ( Cooray fc Kesden|2002 L 

There are many different models of inflation, and they 
often predict very similar power spectra that are close to 
being scale invariant. For models that lead to a measurable 
bispectrum, however, this degeneracy can often be broken 
by studying the dependence of the bispectrum amplitude 
on (£1,^2,^3), e.g. a large signal for squeezed triangles is in- 
dicative of slow-roll multi-field inflation models, equilateral 
triangles are enhanced for models with non-canonical kinetic 
terms, and folded triangles for non-standard vacuum initial 



conditions (e.g. Bartolo et al. (20041, Renaux-Petel (20091, 



Lewis (2011 1). The most studied and constrained form of 



non-Gaussianity is the local ansatz parametrised by a fac- 
tor /nl, and predicted by several inflation models. Current 
constraints on local non-Gaussianity show that the CMB is 
consistent with being Gaussian at the 95% C.L. (Komatsu 



et al.| 2011 1, and constraints on other shapes all show con- 
sistency with Gaussianity as well. 

Given the current limits on primordial NG, astro- 
physical signals are the dominant contribution to non- 
Gaussianity. While Galactic emission and resolved sources 
may be accounted for by masking, unresolved sources and 
residuals have to be characterised. As opposed to primor- 
dial NG, radio sources NG has been detected, and was 
studied for the WMAP mission showing that it yields a 
non-zero flat angular bispectrum parametrised as bsic or 
&DS. The WMA P5 best estimate in th e Q band is &arc = 

Characterising the 



4.3 ± 1.3^it:^sr^(|Komatsu et 



al. 



20091 



NG signal from astrophysical components and more impor- 
tantly from extra-galactic point sources is important for two 
main reasons: (i) to avoid mistaking it for a primordial con- 
tribution (and to allow the development of robust methods 
to isolate primordial NG) and (ii) to learn more about astro- 
physical processes, i.e. to go beyond the description of point 
sources by their number counts and their power spectrum. 
The study of NG from extra-galactic point sources has 
been pioneered by Argueso et al. (20031, focusing mostly 



on the radio sources and including clustering. They showed 
that the point-source angular bispectrum is mostly flat at 
WMAP frequencies and dominates the CMB bispectrum in 
most configurations. [Babich fc Pierpaoli| ( |2008[ ) investigated 
the bias on the /nl estimator induced by radio sources, con- 
sidering the modulation of their number density with ISW, 
of their magnification with lensing and of the fiux cut-off 
with selection effects. Serra fc Cooray ( 2008 1 studied the 



bias on /nl due to radio sources, SZ-lensing and ISW-lensing 
bispectra, investigating the dependence of this bias with the 



resolution scale. Finally Munshi et al. (20091 defined skew- 
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spectra for cross-correlation analysis, derived estimators for 
the skew-spectra and applied it to secondary anisotropies. 

In this paper we study the non-Gaussianity produced 
by infrared and radio point sources in the frequency range 
of the CMB from 30 to 350 GHz, based on numerical simu- 



lations by Sehgal et al. (20101. We investigate the frequency 
and configuration dependence of the angular bispectrum. We 
particularly focus on the non-Gaussianity from IR sources 
and their clustering term. We restrict the study to the sim- 
plest case of full sky maps without masks. Furthermore, we 
do not take into account noise and beam effects. Statistical 
isotropy of all fields considered will be assumed throughout 
this article. The case of noisy, partially masked maps will be 
addressed in future studies. 

In section p] we review the estimator of the (binned) 
angular bispectrum and /nl and develop a parametrisation 
of the bispectrum to display and visualise it efficiently. In 
section [3] we develop a prescription to infer the bispectrum 
from the power spectrum for clustered sources and for differ- 
ent populations. In section |4] we use publicly available full- 
sky simulations of radio and infrared sources to compute and 
characterise their bispectrum at CMB frequencies and we 
compare them to the predictions from the prescription. We 
examine the configuration dependence of the point-source 
bispectra and study the bias they induce on the estimation 
of the primordial local non-Gaussianity in section [5] We fi- 
nally conclude and discuss our results in section [6] 



2 THREE-POINT NG ESTIMATORS 

2.1 Full-sky angular bispectrum estimator 

Given a full-sky map of the temperature fluctuations Ar(n) 
of some signal, it can be decomposed in the spherical har- 
monic basis 



din 



d^n F;„(n) Ar(n) 



(1) 



with the usual orthonormal spherical harmonics Ytm 

/ (fnYimin) Yt*,^,{n) = 5«'(5„,„'. 

Observational data is pixelised, so that the integral is re- 
placed by a sum over pixels. We will assume that the solid 
angle of a pixel, Qpix, is a constant, which is for example 
the case for the HEALPiMjpixelisation scheme that we will 
adopt for the numerical calculations. In this case we have 
that 



(2) 



aim = ^ Ylm{'n^) Ar(nj) ^pix 



This discreteness effect will be important e.g. in section |37L| 
In order to compute the angular bispectrum, which is 
the harmonic transform of the 3-point correlation function, 
we will resort to scale-maps as deflned byjSpergel fc Gold-| 
berg (19991 and also used by Aghanim et al. (20031 and 
DeTroiaet"aI1 ( [2003| ), 



Tt{n) = ^a,„Y«„(n) = f d^n' Ar(n') P,(n ■ n') 



(3) 



^ http://healpix.jpl.nasa.gov 



where Pi is the Legendre polynomial of order £. 



The optimal bispectrum estimator is then ( Spergel & Gold- 



berg|1999 1 



(2^i + l)(2^2-l-l)(24 + l) \0 

d2nr,,(n)r,,(n)r,3(n) (4) 



or it can be written in the form: 



bii 



t2l3 



in 



(2£i + l) (2^2 + 1) (2^3 + 1) 
" h £2 £3 



£2 




E 



mi ,m2 ,^713 



nil 



^•2 

m2 



1713 



^iirni ai2m2 ^^37713 \^ ) 



where the expression in brackets represents the Wigner 3j 
symbols. Equation ([5| is computationally expensive when 
implemented at high £ due to the large number of Wigner 
symbols to calculate. Equation (Ml still requires a few cpu- 
days for a full computation at a Planck-like resolution, 
Nside=1024 - 2048. Binning the multipoles in £, as [Bucher] 
et al. (20101, has the advantage of speeding up the com- 



putations and smoothing out the variations due to cosmic 
variance. 

For a given triangle in harmonic space (£1,^2,^3) the 
number of independent configurations on the sphere is: 



.^ (2^1 + 1) (2^2 + 1) (24 + 1) M £2 
^'-'-'^ = 4^ i 



(6) 



When multipoles are binned in bins of width A£ the 
expression for the scale-maps (Eq. [3| becomes: 



7A^(n) = ^ aimYim{n) 

i£Ai.m 



(7) 



and a binned bispectrum estimator identically weighting tri- 
angles is given by: 



bAei,Ai2At3 ~ 



where 



1 



iVA(A^i,A^2,A4) 



d^nTAii (n) Tai:, (n) TAi, (n) 
(8) 



iVA(A^i,A£2,A4)= Y. ^(^<^2i3 

iiSAii 

One can easily check that the obtained binned bispec- 
trum estimator is unbiased for a constant bispectrum and 
that the bias can be neglected as long as the bispectrum 
does not vary significantly within a bin A£. In the following, 
we have chosen ^max ~ 2048 and a bin width A£ — 64 for 
simplicity and computational speed while retaining enough 
information on the scale dependence (Bucher et al.|2010l. 



2.2 /nl estimator 

The most studied and constrained form of primordial 
non-Gaussianity is the local ansatz, whose amplitude is 
parametrised by a non-linear coupling constant /nl: 



$(X) = t>G{x) + /nl i'i'U^) - {^U^))) 



(9) 



where <I>(2;) is the Bardeen potential and $0(2^) is a Gaussian 
field. This form of NG yields the following CMB angular 
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bispectrum (Komatsu fc Spergel||2001 1 



loc f 

tlt2i3 - 



r dr a^i (r) /3^2 (r) /3f3 (r) + perm. 



(10) 

(11) 
(12) 



with 

2 f 2 
ae{r) = ^ / '^ <ik gT,e{k) ji{kr) 

Hr) = ^ I k\\kP{k)gTAk)Mkr) 

where qt.i is the radiation transfer function, which can be 
computed with a Boltzmann code, je are the spherical Bessel 
functions, and P{k) oc fc"=~* is the primordial power spec- 
trum, with a spectral index Us- 

On large angular scales, the Sachs- Wolfe (SW) effect is the 
dominant contribution to the CMB signal. In this regime, 
the CMB bispectrum takes the following analytical form: 



7 loc 

0^1 £2 «3 °^ 



+ 



p2 «2 ' «2 «2 
<-l •'-2 -"-l <-3 



+ 



p2 02 
'■2 *3 



(13) 



This bispectrum is maximal when one of the multipoles is 
minimal (l\ <^ £2 — £3) which is called the squeezed config- 
uration. 

A commonly used cubic estimator of /nl has been de- 
veloped by Komatsu et al. (20051. It is much faster than 



performing the whole bispectrum analysis and fitting the 
local bispectrum. In its original version, the estimator takes 
into account beam profile and homogeneous noise, and has 
been used on WMAP data to yield the current constraint 
— 10 < /nl < 74 ( Komatsu et al.||2011 1. The estimator was 
then further developed by several authors by adding a lin- 
ear term accounting for masking and inhomogeneous noise 
( Creminelli et al.|2006 1. Here, we will only consider noiseless 
full-sky maps without beam smoothing so that we can apply 
the estimator in its original form. 

To build the /nl estimator we first define the filtered 
maps at comoving distance r and direction n: 



A{r,n) = J2 



aejr) 



a-im y«m(n) 



B{r,n) = -y^at^Ye^in) 



(14) 
(15) 



where Ci is the CMB power spectrum. Bir, n) is then an es- 
timated map of the primordial potential fluctuations $(rn) 
via Wiener filtering. 

The integral of AB^ permits us to estimate /nl as: 



/- E 



B, 



loc \ 



ti^i^^H 



Ct^Ci^Ci^ 



r drd nA(r,n)B^(r,n) (16) 



where B, 



loc 



fj^2«3 \/^i\iit-i ^iii2f--i '^^ ^^^ local bispectrum for 
/nl = 1, to be compared with the observed value Bf^l^^ . 
It can be shown that this estimator takes analytically the 
form: 



/nl 



^tl^l2^t3 Ci Ci Ct 



(17) 



It is near-optimal in the sense that it minimizes the y^ for 
weak NG (under some assumptions on isosceles and equi- 
lateral triangles). The estimator becomes sub-optimal (e.g. 



Eisner & Wandelt (2009 1 ) for large enough /nl, when the 
variance of the bispectrum gets 0(/nl) correction compared 
to the weak NG computation with Wick's theorem (see Ap- 
pendix p|. 



2.3 Parametrisation of the angular bispectrum 

Several ways of visualising the angular bispectrum have been 
proposed in the literature, e.g. isosurfaces in the (^1,^2,^3) 
3D space by Fergusson fc Liguori| ([2010), or slices of constant 
perimeter in the orthogonal transverse coordinate {£±a,£±b) 
space by |Bucher et"aLl ( |2010[ ) . 

Under the assumption of statistical isotropy, the bispec- 
trum htii2t3 is invariant under permutations of £1, £2 and 
^3, ie. it is a function of the shape and size of the triangle 
(^1,^2,^3) only, i.e. independent of its orientation. There- 
fore, we can find a parametrisation invariant under permu- 
tation of ^1, £2, and ^3, which avoids redundancy of informa- 
tion and allows convenient visualisation and interpretation 
of data. Let us denote as (^1,^2,^3) the equivalence class of 
the triplet under permutations. 

The elementary symmetric polynomials ensure the in- 
variance under permutations: 



• (71 = £1 -F £2 + 4 

. (J2 = £i£2 + £i£s + £2£i 
• 0-3 = £1^24 

Through Cardan's fornmla, there is a one-to-one correspon- 
dence between (£1 , £2 , ^3 ) , defined by the roots of the polyno- 
mial X^ — aiX^ + a2X — 0-3, and the triplet (cri, 0-2, (T3). We 
further define the scale- invariant parameters a2 = 12a2/o-l~ 
3 and as = 27a:i/af with coefficient chosen so that a2 and 
5-3 vary in the range [0,1]. As illustrated in the upper panel 
of Fig. [1] this parametrisation does not allow us to discrim- 
inate efficiently between the different triangles. 

We find that the parameters noted (P, F, S) and defined 



• P = cri 

• F = 32(a2 - a3)/3 + 1 

• 5 = 5-3 

provide a clearer distinction of the triangles as is illustrated 
in the bottom panel of Fig. [T] 

To illustrate the use of our (P, F, 5')-parametrisation, we 
plot in Fig. IWlthe theoretical CMB bispectrum produced by 
the local NG model /nl, computed through Eq.flO] The tri- 
angle perimeters, P, vary between Pmin = 30 (equilateral 
configuration with Imin = 10) and Pmax = 6120 (equilat- 
eral configuration with ^max = 2040). We plot representa- 
tive perimeters tracing the building up of the bispectrum 
with scale, giving P/3 on each panel. The color code scales 
from deep purple/black (most negative) to red/dark grey 
(positive). 

In the first panels for the smallest perimeters, few trian- 
gles are present. As the perimeter increases the (P, S) space 
is populated starting with equilateral configurations (upper 
right corner) first to reach squeezed configuration (upper left 



^ For space reasons we only include some of the bins in the figures 
in the paper. Full resolution plots with all perimeter bins are 
available upon request. 
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Figure 1. Triangles of constant perimeter, P, in tlie parametri- 
sation defined by the normalised symmetric polynomials (upper 
panel), or in the parameter space (F, S) defined in the text (bot- 
tom panel). 




Figure 2. Theoretical CMB angular bispectrum produced by the 
local model in the (P, F, 5)-parametrisation. With chosen repre- 
sentative perimeters and asinh color scale from deep purple/black 
(most negative) to red/dark grey (positive). P/3 is given on each 
panel. The axes are the same as those in Fig. ^ 



corner) later. Conversely in the last panels for the largest 
perimeters, the resolution limit ^max = 2048 limits the pos- 
sible configurations, leaving only the equilateral triangles in 
the end. 



The Sachs- Wolfe shape (see Eq. 131 is visible at low perime- 
ters, with the colors (value of the bispectrum) varying hori- 
zontally with S but not vertically with F. We note that the 
strong negative values (deep purple/black) are located in the 
near-squeezed configuration (upper left corner). The sign of 
the radiation transfer function can be traced via the equilat- 
eral triangles which are positive for P/3 ~ 200 (first acoustic 
peak) and become negative for P/3 ~ 500 (second acoustic 
peak) etc. The squeezed configuration in turn is always neg- 
ative as the smallest multipole has negative transfer function 
via SW. The folded configuration has relatively large nega- 
tive values when the two smallest multipoles are in the first 
acoustic peak while the biggest multipole is in the second 
acoustic peak (P ~ 900 - 1200 => P/3 ~ 300 - 400). For 
larger perimeters the structure becomes complex as several 
acoustic peaks intervene. 



3 PHYSICAL PRESCRIPTION FOR THE 

POINT-SOURCE ANGULAR BISPECTRUM 

In this section we develop a prescription which permits us 
to predict the (bi)spectrum of point sources, starting from 
the simplest case of a single randomly-distributed popula- 
tion to the case of multiple clustered populations. This is 
a situation that is encountered in current and future CMB 
analyses. Indeed until recently, CMB experiments have fo- 
cussed on frequencies where unclustered radio sources are 
the only dominant kind of point sources, but the CMB is 
also non-negligible at higher frequencies where an indepen- 
dent population of dusty galaxies becomes important to- 
gether with the SZ signal of clusters. This is of particular 
relevance for Planck which has a large frequency range cov- 
ering both populations. 



3.1 Single source population: shot-noise 
contribution 

A source with flux S enclosed in a pixel with solid angle 



pix yields a temperature variation AT 

aB(v,T) I 



k^T^—, where 



kv = — Qj T B{u, T) is the black-body spectrum and 

TcMB is the CMB mean temperature. 

As shown in AppendixlAlthe power spectrum of a source 
population is given by: 



Ct 



^clust I ^shot 



(18) 



The discreteness of the sources produces a constant-term 
spectrum C*°* which is usually named 'Poissonian' because 
the number of unclustered point sources is driven a priori by 
Poisson statistics ( Sehgal et al.|[20T0 |. The shot-noise term 
reads: 






47r 



9'— d9 



(19) 



where ^ is the differential number counts of sources and 
Scut is the detection limit, i.e sources with S > Scut are 
detected and masked, the rest being unresolved. 
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The discreteness property of sources, when computing 
the three-point correlation function, yields a statistically 
isotropic angular bispectrum constant with £: 



TShot 



An 



9'— d9 



(20) 



for £i / 0. Equations H8f, u9\ and pOl) are derived in more 
detail in Appendix [X] 

The case of sources randomly and independently dis- 
tributed on the sky is that of the radio sources at 
CMB frequencies. The correlation vanishes and the total 
(bi)spectrum is equal to the shot-noise (bi)spectrum. The 
distribution of the sources is that of a white-noise entirely 
characterised by the one-point probability distribution. In 
this case, Eqs. (191 and (201 for the shot noise contribution 



can be reformulated simply in terms of temperature vari- 
ations where the white-noise spectrum and bispectrum are 
related to the variance a^ and skewness Kg of AT: 



C white 2 o 
£ — a iipix 



6Avhite (-x2 

^1^2^3 — ^3 "pi; 



(21) 
(22) 



with a 



{{AT - {AT)f) and K3 = {(AT - {AT)f). 



3.2 Single source population with correlations 

The effect of clustering of a single population of sources, 
namely radio sources, on the bispectrum was pioneered by 



Argueso et al. ( 2003 1 who proposed a prescription to address 



this issue. The elements entering the prescription are the 
number counts of sources and a theoretical or observational 
correlation function w{6). Defining 



(23) 



P(fc)ciust = 27r / w{e)Jo{ke)de 



where 9 is the distance on the flat patch and Jo the Bessel 
function of the first kind and of order zero, the prescription 
is: 



where (5g are the Fourier coefficient of the map, 'shot' and 
'white' refer respectively to the shot-noise and white-noise 
process. Then Argueso et al. compute the angular bispec- 
trum from simulated square maps based on their prescrip- 
tion. 

In our study, we have extended the above-described pre- 
scription to analytically derive the full-sky angular bispec- 
trum. The full-sky scale-maps read: 



rr(n) = 






CT 



xt; 



= (n) 



which reduces to Eq. ( 24 1 in the flat-sky approximation £ — 
2nk, and assuming P{k) does not vary much within a k bin. 
So the power spectrum is given by: 



C, 



tot y'-^clust I ^sliot 



We remind the reader that the Argueso et al.'s prescription 
is speciflcally made to meet the above relation, and that 



More interestingly the bispectrum is found to be: 



"«l<?2f3 



Ci 



/^ white 



1 + 



/^clust 

/^whitc 
^12 



CI 



/"fwhite 



X bp 



(25) 



/-ftot /-<tot /^tot 
^f, ^lo '-■I:, 



with bj^i^i^ = const and a = 



/^whii 



Equation (251 of the bispectrum relates to the clustered 



power spectrum, or conversely the correlation function en- 
tering in the prescription. 

3.3 Two populations of sources with clustering 



The previous prescription, Eq. (24 1, describes well a single 



point source population. However, it fails at describing the 
cases where more than one population contribute to the sig- 
nal. In particular, the case of the present generation of CMB 
experiments observing the CMB from low (30 GHz) to high 
(860 GHz) frequencies calls for an appropriate prescription 
which deals with independent source populations. Indeed, 
early results from the Planck mission ( Planck-Collaboration 



2011c I (Planck-Collaboration 2011a I show that radio and 



IR galaxies both contribute at frequencies of 100 GHz and 
above. In the following we thus extend and generalise the 
prescription accordingly. 

We hence define the harmonic coefficients as: 



(tot) 



fwliitc.l) , 

aL + 



/-*clust _] /^shot 

^i + *-"<; (wiiitc,2) 

^white,2 em 



(26) 



1 (white, 1) 1 (white. 2) .1 n . T , • 

where a^^ and a^^ are independent realisations 

of white-noise with different number counts. Index 1 refers 
to the radio population and 2 to the infrared population. 
The spectrum has then a form similar to the case of a single 
source population: 

^ ^^-fclust I ^.-*shot,l+2 

But the angular bispectrum differs, and reads: 



01112(3 



T white,! 1^ 

' °iii2e3 + 



\ C7i 



white, 2 



c; 



c: 



^white,2 f~i 



3 y r white, i 

white, 2 ^ "^l<?2<3 



(27) 

For illustration, let us compare equations ( |25[ ) and ( |27[ ) 
in the equilateral configuration with white-noises of both 
populations derived from the same number counts, and ne- 
glecting the shot-noise of the second population: 

3/2 



Llpop 



c? 



/-fwliitc 



,2pop 



1 + 



c: 



clust \ 3/2 



/^ white 



(28) 



(29) 



The difference between these two formula is maximal when 



clust //^^vhite 



C clust //-I- 
e /^i 



1 and can be up to 40%, illustrating the 



need to properly account for the different populations. 

The two-population case, representative of the CMB 
context in the frequency range of interest, can be straight- 
forwardly generalised to more point-source populations with 
or without clustering properties. 
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4 RESULTS 

In this section we present the bispectra of radio and IR 
sources computed on simulations described below, showing 
the configuration dependence of the angular bispectrum and 
its frequency behaviour. 



4.1 Data used 

For our analysis, we used the all-sky simulated mapaj of 
the IR and radio point sources provided by [Sehgal et al.| 
(20101) at 30, 90, 148, 219 and 350 GHz. We provide here 



Equilateral 



Orthogona 



a brief summary of the map construction. For a detailed 
description, we refer the reader to the above-cited article. 
The maps are based on N-body simulations of the large scale 
structure, with a volume 1000 h~^Mpc on a side, produced 
using a tree-particle mesh code. Dark matter (DM) haloes 
are identified and are then populated with infrared and radio 
galaxies. The model for the radio sources is adapted so that 
the radio luminosity function matches that of the observed 
radio sources at 151 MHz. 

The IR-source model was partially based on|Righi et al.| 



( 2008 1. The DM haloes are populated with galaxies of given 



luminosities taking into account a Poisson term and a corre- 
lation term. The model was constructed so that it is compat- 
ible with constraints on the luminosity function, the source 
counts and the fluctuations from SCUBA, BLAST, Spitzer 



and ACBAR (see Sehgal et al. (20101 for details). However, 



the simulation of IR sources, used here, does not account for 
the most recent observational constraints from ACT, SPT 
Herschel and Planck results. 

Maps of the different astrophysical contributions, in 
HEALPix format at Nside=8192, were produced by repli- 
cation of one simulated octant of the sky. This procedure 
does not properly account for the signal at the largest scales 
especially up to the octopole, ^ = 3 ; it also introduces ex- 
cess power at ^ ^ 300 for infrared sources as discussed by 



Sehgal et al. (20101, but the power is correct for higher mul- 



tipoles. The produced maps have half-arcminute resolution 
but we degraded them to Nside=1024 and used a uniform 
binning A£ = 64, to speed up computations. We checked 
that this procedure does not introduce a bias by comparing 
the binned spectrum in the degraded map to the unbinned 
spectrum in the original map and flnd excellent agreement. 
The octant replication in the map- making translates mainly 
into a lack of power in the first bin (centered around £ = 32) 
which is hence discarded in later plots. 



4.2 Radio source characterisation 

We first investigate the bispectrum dependence on the con- 
figurations at a given frequency. We plot the angular bispec- 
trum in four commonly used configurations, namely equi- 
lateral {£,£,£), isosceles orthogonal {£,£,\/2£), isosceles fiat 
{£,£,2£), and squeezed {iinin,£,£) configurations. They fur- 
thermore sample rather well the configuration space (see Fig. 



0. 



We see in Fig. |3J black thick line, that the bispectrum is 



^ The frequency maps are available on Lambda website. 




100 1000 

multipole I 



100 1000 

multipole 



Squeezed 




http://lambda.gsfc.nasa.gov/tooIbox/tb_cmbsim_ov.cfm 



Figure 3. Thick black line: Radio angular bispectrum at 90 GHz 
in different configurations: equilateral {£,£,£), isocele orthogonal 
(£, £, V2£), isocele flat {£, £, 2£), and squeezed (£min, ^, £)■ The thin 
red line is the fit with prescription. The unit of the bispectra is 
Jy'^ ■ sr~^ 



constant. This result is independent of the frequency. More- 
over the value of the constant is independent of the configu- 
ration. This is what we expect from white noise and it shows 
that radio sources are indeed randomly distributed over the 
sky. 

We show in Fig.|4]the dependence with frequency of the 
bispectrum amplitude averaged over all the configurations, 
in equivalent temperature unit for the upper panel and fiux 
unit for the bottom panel. The amplitude is maximal at the 
lowest frequency 30 GHz, it then decreases to become mostly 
fiat above 90 GHz because of free-free emission and inverted 
spectra sources. 

The bispectrum was fitted with the one population pre- 
scription described in Sect. |3.2[ using the power spectrum 
extracted from the simulations and the multiplicative con- 
stant which minimizes the x^ of the observed bispectrum to 
the prescribed one with error bars from Wick's expansion 
(see Appendix [B| . The fit of the bispectrum with the pre- 
scription is very good: Fig. |5] shows that the relative error 
(exact to fit) lies in the range between -2% and -1-2% with a 
mean relative error always less than l%o. 

4.3 IR source characterisation 

Figure [6] shows the amplitude of the angular bispectrum 
in the equilateral configuration for the different frequencies. 
Dusty galaxy emission peaks at IR frequencies, and plum- 
mets in the radio domain and so does the amplitude of the 
bispectrum. The bispectrum decreases with £, well fitted by 
a power-law for 100 ^ £ ^ 1000, and with a fiattening at 
higher multipoles. The decrease is expected from the clus- 
tering of IR galaxies on large scales. We found that the slopes 
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Figure 4. Amplitude of tlic radio angular bispectrum, averaged 
over all configurations, at 30, 219, 148, 277, 90, and 350 GHz, 
in temperature units for the upper panel and fiux units for the 
bottom panel. 



were about the same for equilateral, flat and orthogonal con- 
figurations except for the squeezed triangles. Indeed, for the 
latter £i is fixed, while for the other configurations £i,£2,f.3 
are all proportional to £. 

The flattening of the bispectrum at high multipoles, in- 
dicative of the shot-noise contribution, occurs at lower mul- 
tipoles with increasing frequency (e.g. at 350 GHz the bis- 
pectrum deviates from a power law at ^ ~ 1000 while at 
90 GHz it is not before £ ~ 1500). This is explained by the 
contribution of the high-flux galaxies, in Sehgal et al.'s simu- 
lations, that accounts for the shot-noise [b'''^°*- ~ J S^f^dS) 
and at the same time have a steeper emission than the galax- 
ies accounting for the clustering term. 

The bispectrum was fitted with the one population pre- 
scription described in Sect. |3.2[ using the multiplicative 
constant which minimizes the x^) ^-s previously for radio 
sources. We show in Fig. [T] how the bispectrum of the IR 
sources compares with the prescribed one. We see that the 
bispectrum obtained with the prescription is good, with a 
mean relative error always ^ 5% and a standard deviation 
^ 31%. At 350 GHz an outlier at -400o- was discarded for 
the computation of the standard deviation. Figures [5] and IT] 
show the distribution of the relative error between the bis- 
pectrum derived from the prescription and the actual bis- 
pectrum measured in the simulated maps. In other words, 
it exhibits departures from the predicted bispectrum values. 
The dispersion of these relative errors is larger for the IR 
sources (Fig. n\ than for the radio sources (Fig. [5| at all the 
frequencies. This behaviour is not an indication of a mis- 
match between the predicted and the actual bispectra; it 
relates to the intrinsic sample variance of both the IR and 
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Figure 5. Relative error distribution with the prescription for 
radio sources alone. The vertical bar shows the mean of the dis- 
tribution 
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Figure 6. Infrared equilateral bispectrum at 90, 148, 219, 277, 
and 350 GHz from bottom to top. Dotted corresponding lines are 
the power-law fit 
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infrared sources alone. The vertical bar shows the mean of the 
distribution. The quoted standard deviation at 350 GHz is com- 
puted after discarding one unique negative bispectrum outlier at 
400cr. 



radio bispectra. The IR sources bein g weakly non-Gaussian 
(the value for a defined in Eq. ( 25 1 is 3 x 10"^ for the IR 
sources, compared to 0.3 for radio sources) the variance of 
their bispectrum is indeed large compared to the bispectrum 
value (see Appendix [b|) . 



4.4 Total contribution from IR and radio source 
populations 

We now present the results when the two populations of 
sources contribute to the signal at the frequencies 30, 90, 
148, 219, 277, and 350 GHz. To do so, we simply add the 
simulated maps at each frequency. 

We illustrate the angular bispectrum dependence on fre- 
quency for one single configuration, namely equilateral, see 
Fig-i 

The frequency behaviour is as expected from an inde- 
pendent combination of the IR and RAD bispectra. The 
radio source contribution dominates at low frequencies 30 
and 90 GHz (blue and purple lines) and its bispectrum is 
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Figure 8. Infrared -|- radio orthogonal bispectrum at 30, 90, 148, 
219, 277, and 350 GHz. The purple fiat upper line is 30 GHz, the 
blue flat lowest line is 90 GHz, more variable decreasing lines from 
bottom to top are 148, 219, 277, and 350 GHz 



fiat. Infrared galaxies dominate at the highest frequencies 
277 and 350 GHz (black and red upper lines) and show the 
characteristic power-law dependence due to clustering fol- 
lowed by a fiattening of the bispectrum. At intermediate 
frequencies both populations contribute to the signal. The 
clustering-induced term of IR-galaxies dominates on large 
angular scale while the random-noise term of radio-galaxies 
dominates at small angular scale. The cross-over between 
radio and IR-galaxy bispectra is shifted to higher Is with 
increasing frequency. It is worth noting in Fig. |9] that at the 
lowest multipoles and at highest frequencies, the IR galaxies 
produce a bispectrum at least 10 times more important than 
the radio sources. 

We illustrate the angular bispectrum dependence on 
configurations at 350 GHz, see Fig. [lOlthick black line. The 
errors bars were computed in the weak NG approximation 
see Appendix [B] We note that infrared-radio cross-over oc- 
curs at about the same scale for the equilateral, orthogonal 
and fiat configurations, but is at higher i in squeezed con- 
figuration. This is expected because the squeezed IR bispec- 
trum decreases more slowly than other configurations, one 



of the multipoles being fixed. Figure 10 also displays (thin 



red line) the bispectrum computed with two-population pre- 
scription derived in Sect. [3] i.e. adding up independently the 
prescription for radio sources and infrared sources derived 
in the previous sections. This is compared to a bispectrum 
computation considering only a single population (thin blue 
line). From Fig. 10 and Fig. 11 showing the distribution 



and Fig. 

of relative errors with respect to these two prescriptions, 
we see that the two-population prescription performs much 
better than the single-population prescription. The former 
captures well the overall shape of the bispectrum and it ad- 
justs particularly well the high ^s. As a matter of fact, the 
mean relative error is lower than 2.5% up to 350 GHz and 
the dispersion increases from 1% at 30 GHz to 21% at 350 
GHz. 

As expected, at 30 and 90 GHz the two prescriptions 
give same results since radio sources totally dominate the 
signal. At higher frequencies, both the mean errors and the 
dispersions derived using the two-population prescription 
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Figure 9. Ratio of tiie IR angular bispectrum to tiie radio one at 
90, 148, 219, 277, and 350 GHz from bottom to top. The dotted 
line indicates equality. 



are smaller than those obtained with the single-population 
prescription. Interestingly enough, at the highest frequency 
(350 GHz) where IR emission from galaxies is dominant, 
the single-population prescription is not satisfactory. As a 
matter of fact, configurations with at least one high multi- 
pole dominate the distributions (e.g. 7/8th of the configu- 
rations have at least one ii ^ 1000). At 350 GHz these £s 
are dominated by infrared shot-noise, so the computation of 
the prescriptions combine the IR shot noise and radio spec- 
trum which are both flat. The radio emission is subdomi- 
nant compared to infrared shot-noise but nevertheless not 
negligible so the single-population prescription leads to an 
overestimate of the total bispectrum. This is clearly visible 



in Fig. 10 where the single population prescription (blue thin 



line) is systematically higher than the computed bispectrum 
(black thick line) and than the two-population prescription 
(red thin line), particularly at high multipoles. 



5 CONSEQUENCES ON NON-GAUSSIANITY 
MEASURES 

5.1 {P,F,S) parametrisation 

The parametrisation of the angular bispectrum proposed in 
Sect. |2.3| allows us to visualise the bispectrum dependence 
on the configurations. The bispectrum is computed for 32 
perimeters in the I space. Only nine perimeters are shown 
for illustration in Figs. [12] and |13| for the radio and IR- 
source populations respectively. The bispectrum values are 
colour coded from blue (lowest value) to red (highest value). 
The succession of plots, arranged by increasing perimeters, 
exhibits the allowed configurations at given perimeter with 
the equilateral configuration being the starting (upper left 
panel) and ending point (lower right panel). Unsurprisingly 



Figure 10. Infrared -|- radio bispectrum at 350 GHz in differ- 
ent configurations, with error bars. The thick black line is the 
computed bispectrum. Smoother solid lines are the fit with pre- 
scriptions: the upper blue line with the single-source prescription, 
and the lower red line with the two-sources prescription. 



the bispectrum amplitude of the radio sources does not vary 
with the configuration (same color code for all points in Fig. 
121. ^ 



13 the bispectrum ampli- 



As for the IR sources. Fig. 
tude decreases with perimeter, thus from upper left to lower 
right panel. Moreover, it is worth noting that the bispec- 
trum values do not vary vertically. This means that within 
the proposed (P, F, S')-parameterisation the bispectrum of 
the IR sources is quasi-independent of F, reducing the full 
bispectrum to a function of the two parameters P and 5*. 
Finally, at a given perimeter, i.e. scale, the bispectrum of 
the IR sources is more dependent on the configuration and 
peaks in the squeezed triangles, upper left points in panels 
2 to 5. 



5.2 point sources contamination of /nl 

We now explore the point-source non-Gaussianity in terms 
of contamination of the /nl estimation. For pedagogical pur- 
poses, we consider the SW regime, i.e. a constant transfer 
function without acoustic oscillations, damping etc. 

For IR sources alone, which dominate at high frequen- 
cies, 6(^1,^2,4) oc yCViTT^^C^ 
latest constraints on the GIB 



Gollaboration|201 Id I , this yields: 

2, £3) OC 



Gombined with Planck's 
I X Ce ~ const (Planck- 



b'^'iiu 



VhMi 



which has a similar shape to the local template in the SW 
limit, Eq. (13 1. We define the scalar product between the 



bispectra /, g: 
<f,9>= 






f(£l,l2j3)xg{h,i2j3) 
/T-iClVIB /^CMB r<CMB 

^tl ^t2 ^ts 
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Figure 11. Relative error distribution with the prescriptions for 
the combination of infrared and radio sources. The dotted line is 
obtained with the single population prescription, the solid line 
with the two populations prescription. Vertical bars show the 
mean of the distribution, blue for the single population prescrip- 
tion and red for the two populations prescription. 
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Figure 12. Radio bispectrum in the (P, _F, S) parametrisation at 
90 GHz. Each plot is a slice of constant perimeter P, the value 
of the bispectrum is encoded in a logarithmic color scale from 
violet-blue to red. The axes are the same as those in Fig. [l] 
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Table 1. Bias on the /nl estimator, A/^*^ , due to radio sources 



where the denominator is the variance of the local bispec- 
trum, for triangles with (.\ ^ I2 i^ (■?,■ The correlation coeffi- 
cient between a bispectrum 6° and the local bispectrum 6'°'^ 
is: 



< 6", 6'°'= > 



cos Va — 



w° 



with a being IR or radio. 

Figure [14] shows that, for £ < 200, the correlation be- 
tween radio and local bispectra decreases quickly, so that the 
two bispectra may be distinguished efficiently. Conversely, 
even when using a large multipole range, the IR bispectrum 
is significantly correlated with the local one. 

The contribution of the bispectrum of a point source 
population a to /nl is: 



A/5l 



< fe",6^«L= 



Il6"li 



(30) 



< 5^NL=i^ 6/nl=i > ||&/nl=1|| 
This equation is the usual bias ( Serra fc Cooray|2008 1 of the 



local-optimised NG estimator, when the local bispectrum 
has the form of Eq. ( 13 1. 



A more comprehensive computation of A/^l is achiev- 
able by applying the full local estimator described in Sect. 
2.2| We built up this estimator using the full transfer func- 



tion from the latest version of CAMB ( Lewis et al. 2000 1 



with WMAP7-I-BAO-I-H0 cosmological parameters (Larson 
et al.|20rT|), and we tested the estimator on simulations from 



Eisner & Wandelt ( 2009 \ . We found the previously noted re- 



sult that the variance of the estimator increases with /nl. 
It is unbiased in the range we have tested (0 ^ /nl ^ 200). 
We used this /nl estimator on two sets of simulated 
maps: maps containing all the point sources, and maps with 
only sources below the flux limit of Planck's Early Release 
Compact Source Catalogue (ERCSC) ( |Planck-Collaboration| 
2011bl), namely 5c =0.5, 0.5, 0.3, 0.3, 0.3, 0.25 Jy as a func- 



tion of frequency. Moreover, we have computed the estimator 



at three resolutions, 



recalibrating the Sprim normali- 



sation in each case. Tables [T] and [2] summarise these results. 
The bias A/nl ° (see Table Ilk is negative on large an- 
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Figure 13. IR bispcctrum in the (P, F, S) parametrisation at 148 
GHz. The axes are the same as those in Fig. IT] 




Figure 14. Correlation between SW local bispectrum and re- 
spectively IR (red) and radio (blue) bispectrum, as a function of 
the maximum multipole used (^niin = 2) 



without flux cut 


u (GHz) 


30 


90 148 


219 


277 350 


-t-max ^ OU 
^max = 700 
^max = 2048 


-3.5-10- 
-1.3-10- 
-1.8-10- 


-**-5.9-10-«-9.2-10-s 
-^ -0.00019 -0.0033 
-5 -0.0026 -0.039 


-0.0027-0.023-1.0 
-0.063 -0.55-9.0 
-0.68 -4.8 -67 


with flux cut 


u (GHz) 


30 


90 148 


219 


277 350 


^max = 700 
^max = 2048 


-1.3-10- 
-1.8-10- 


-f" -0.00019 -0.0033 
-5 -0.0026 -0.039 


-0.078 
-0.67 


-0.74 -11 
-6.3 -66 



Table 2. Bias of the /jml estimator, A/^^, due to IR sources 



gular scales, for ^max = 50. The bias due to radio sources 
becomes positive at higher multipoles in agreement with 



Serra fc Cooray| ( |2008| ). This is due to the CMB bispec- 
trum being negative in the SW-dominated regime and the 
radio bispectrum being positive. The bias increases by 5 or- 
ders of magnitude at the highest resolution, ^max = 2048. 
The reason for the rapid increase of the bias with i'max re- 



lates to the weight of the observed bispectrum in Eq. (171, 
B^i'^i^i^/Ci^Ci^Ci^, which rapidly increases with multipole 
as the product of spectra decreases more quickly than the 
bispectrum. This leads to a 1/Ci dependence in squeezed 
configurations and to a 1/C| dependence in equilateral con- 
figurations. When the observed bispectrum is associated to 
CMB signal alone, its decrease cancels the increase of the 



weights so that the sum in Eq. (171 converges. Conversely, 
the sum diverges when the observed bispectrum is associ- 
ated with a non-CMB signal and does not decrease with I 
as fast as the CMB. 

The bias Af^j^ ° is maximal at 30 GHz and rapidly 
decreases with frequency. It slightly increases again at the 
two highest frequencies following the amplitude of the bis- 
pectrum in temperature units which is plotted in the upper 
panel of Fig. ^ The relative error of A/^^° for £niax = 700 
is of the order of 1.5% independently of the frequency. It 
amounts to 2.3% for ^max ~ 2048. These errors bars were 
computed with simulations using the catalog of sources 
present in Sehgal et al.'s maps 

As shown in Tablell] masking sources above the ERCSC 
flux limit proves very efficient to significantly decrease the 
radio contamination to /nl at all the frequencies. At a 
Planck-like resolution, ^max = 2048, the bias A/n]^° is re- 
duced below unity above 150 GHz. It is of the order of 
Planck's expected error bars at 90 GHz. At 30 GHz the 
bias is still important. 

The bias due to IR sources Af^ is always negative, see 
Table [2] As a matter of fact, we have shown that the IR bis- 
pectrum peaks in squeezed configurations just like the CMB 
bispectrum and these configurations thus dominate the sum 



in Eq. (17 1. Moreover, in the squeezed limit the CMB bis- 



pectrum is negative while the IR bispectrum is positive. For 
the same reason as for radio sources, the bias A/^^ blows 
up at high multipoles. This is particularly important at a 
Planck-like resolution, £max = 2048, where primordial NG 
tests will need to carefully handle the contamination by IR 
sources. The IR sources emission plummets at radio frequen- 
cies so that A/n^ is completely negligible below 220 GHz. 
It becomes of the order of Planck's error bars at 277 GHz 
and it reaches WMAP's central values for /nl at 350 GHz. 
The relative error of Af^ ranges between 6 and 7% from 
148 to 350 GHz for ^^ax = 700. It ranges between 3 and 7% 
for .^max = 2048. (These error bars were computed analyti- 
cally with the weak NG approximation - see Appendix [B|) 
At higher frequencies the IR contamination to the bispec- 
trum is likely larger but the contamination from our Galaxy 
needs to be taken into account as well. 

Interestingly, masking sources above the ERCSC flux 
limit does not diminish Af^, as most of the IR sources 
are unresolved and the IR clustering is mostly due to faint 
sources. Masking may even artiflcially boost Af^, for ex- 
ample at 277 GHz, since it mostly affects the flat shot-noise 
which produces a positive bias A/nl. 



Characterisation of the non-Gaussianity of radio and IR point sources at CMB frequencies 13 



6 CONCLUSIONS AND DISCUSSION 

We have studied the non-Gaussianity produced by point 
sources in the frequency range of the CMB from 30 to 350 
GHz. We have developed a simple and accurate prescription 
to infer the angular bispectrum from the power spectrum of 
point sources, considering different independent populations 
of sources, with or without clustering. 

Using publicly available all-sky simulations of radio and 
IR sources, we have computed the full-sky binned bispectra 
for these two populations of sources. We have compared the 
measured bispectra to those predicted from our prescription 
and found a very good agreement between the two. We have 
displayed the angular bispectrum using a new parametri- 
sation which highlights efficiently the configuration depen- 
dence. 

We have characterised the angular bispectrum of the IR 
and radio sources showing the configuration dependence and 
the frequency behaviour. In particular and for the first time, 
we showed that the IR bispectrum peaks in the squeezed 
triangles and that the clustering of IR sources enhances the 
bispectrum values by several orders of magnitude on large 
angular scales £ ~ 100. The bispectrum of IR sources starts 
to dominate that of radio sources on large angular scales at 
150 GHz, and it dominates the whole multipole range at 350 
GHz. 

Finally to illustrate the contamination of local CMB 
non-Gaussianity by point sources, we derive the bias on /nl 
induced by radio and IR sources, for WMAP or Planck-like 
angular resolutions. Radio sources produce a positive bias 
which is significantly reduced (A/nl < 1 for i/ ^ 150 GHz) 
by masking the sources above a given fiux limit taken as 
the ERCSC cut. The form of the IR bispectrum mimics a 
primordial 'local' bispectrum /nl on large angular scales. 
The IR sources produce a negative bias which becomes im- 
portant for Planck-like resolution and at high frequencies 

(A/nl ~ -6 at 277 GHz and A/nl 60-70 at 350 GHz). 

Most of the signal is associated with the clustering of faint IR 
sources. Therefore, the bias A/n^ is not reduced by masking 
sources above a fiux limit but, in some cases, even increased 
due to the reduction of the shot-noise term. 

Our analysis highlights the sensitivity of the bias on /nl 
to the experimental properties (maximum resolution and fre- 
quency range), the point-source models (clustering or no, 
resolved or not) and their scale dependence with respect to 
the CMB. For high resolution high frequency CMB experi- 
ments, primordial NG estimations need to take special care 
of astrophysical contaminations. One solution would be to 
estimate the primordial and astrophysical non-Gaussianity 
simultaneously. 
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APPENDIX A: CLUSTERED SOURCES 
SHOT-NOISE 

A source with fiux S enclosed in a pixel fipix yields a rise of 
temperature compared to the CMB: 



AT: 



S 



(e- 



2!/2fe£ 



(Al) 



BBjv.T) I 
ST 



'Tcr, 



where x — hu/kBTcMB and ki. 

The two-point correlation function of point sources 
takes the form: 

(Ar(n)AT(n')) = F(n, n') + r5„,„, 

where F{n, n') is the correlation function coming from the 
spatial distribution of the sources, and the Kronecker term 
comes from the discreteness of the sources: 

r = {Ar^)-F(n,n) 

Assuming statistical isotropy, we get: 



Ci = cf ""' + cf°' with cf°' = r fipi. 



Indeed: 



(aimae,^,) = Ad^nd^n' y,„(n) y;„,(n'){Ar(n)Ar(n')) 



= CSl"='<5„, 4\, 






Of Offi Oram 



= Ci "°' 5u' 5^ra' + Ffipix / d rvYim{n) Y^*,„/(n) 

= C^ 5ui Sram' + P fJpix S^i 5rnm' (A2) 

Ci has units ^,K^ ■ sr. Let us number by i=l..N all sources 
of the sky, then the temperature of a pixel is given by: 



rCi/ 



Ar(n) = -^^S, xl[,. 



(A3) 



where l[i6n] is 1 if the source i is in the pixel and otherwise. 
We have: 



http://lambda.gsfc.nasa.gov/toolbox/tb_cmbsim_ov.cfm 



{l[ieni) = 1/npix 
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Hence for n 7^ n': 

F(n,n') = (AT(n)Ar(n')) 



When the field is close to Gaussian, the 6-point correla- 



Then we find: 



„-^(^S'iSj X l[,gn]l[jGn']) 



Q^(Z2 ^'^^ ^ l[»6n]l[jen]) 
P'" ij = l 

pi" i^j 



U2 N 



+ r^E^'x<l[»en,) 



Recalling fipix 
of sources 

(ATinf) = 
And finally: 



lim/Ar(n)Ar(n')) + ^— J2 S' 

l-"^^ ^ sources 

^^ and introducinff 4? the number counts 



F(n,n) + 



A-rrQn 






cf- = rn,..= 'l^rs^^AS 



(A4) 



which is the shot-noise formula Eq. ( 19 \ 



The integral runs from S=0 to Scut the fiux detection limit 
of the survey, ie sources with 5 > 5cut have been removed. 
Note that this result is independent of the two-point corre- 
lation function, which we did not specify. 
At order 3 for the angular bispectrum, the computation is a 
bit more involved but follows the same line, and we find: 



{ae^miae^m2ae3m3) = G^f^g™^ x ''<i<2<'3 (A5) 



with 



7 shot 



47r 



dn 
dS 



dS 



(A6) 



which is the shot-noise formula Eq. ( 20 \ 



APPENDIX B: BISPECTRUM VARIANCE IN 
THE WEAK NG APPROXIMATION 



The bispectrum estimator Eq. Q can be put in the form: 

ri23aia2a3 (Bl) 



&123 = j^— X ^ Gl 
JV123 ^ — ' 

mi2.3 



where a shortened notation is used: A''i23 is the number of 
triangles defined in Eq.16] and 6*123 is the Gaunt coefficient: 

G123 = / d^nY«imi(n)Ff2m2(n)yf3„3(n) . (B2) 

;he form 

/ ^ G123 Gi'2'3' 



Then the bispectrum covariance takes the form 

Cov(6i23 ,6l'2'3') = 



^123 ,"^123 

{{aia2a3aiia2'a3,) — {aia2a3){ai/a2'a3/))(B3) 



tion function can be computed with Wick's theorem (Wick 



(19501, Komatsu et al. (20021 and references therein 



{aia2a3ai/a2'a3') = (Gf)i23 (5i*,i' (52*,2' <53*,3' + 14 perm. 

(B4) 
where 5,.,^ = {-l)"^^5i^^i.S-,n,,mj- 

The 15 permutations of (1,2,3,1', 2', 3') are listed below along 
with their contribution 5Cov to the covariance. 



(1*2)(3*1')(2'*3' 
(1*2)(3*2')(1'*3' 
(1*2)(3*3')(1'*2' 
(1*3)(2*1')(2'*3' 
(1*3)(2*2')(1'*3' 
(l*3)(2*3')(l'*2' 
(1*1')(2*3)(2'*3' 
(1*1')(2'2')(3*3' 
(1*1')(2*3')(3*2' 
(1*2')(2*3)(1'*3' 
(r2')(2'l')(3*3' 
(1*2')(2*3')(3*1' 
(1*3')(2*3)(1'*2' 
(1*3')(2*1')(3*2' 
(r3')(2*2')(3*l' 



5Cov = except if £3 = I'l = (B5) 

5Cov = except if ^3 = 4 = (B6) 

5Cov = except if ^3 = ^3 = (B7) 

5Cov = except if £2 = £[ ^ (B8) 

5Cov = except if ^2 = ^2 = (B9) 

5Cov = except if ^2 = ^3 = O(BIO) 

5Cov = except if £1 = l'^ = O(Bll) 

5Cov = 5^i|^i^5,^,, 5,,,, 5,3,, (B12) 



iVl23 

5Cov = ^^4^2^ S^ 

"123 



SCoY = except if ^1=4= 0(B14) 
5Cov = £^i^£3^5,,,, 5,^,, 5,3,, (B15) 



5Cov = 



iVl23 "•^1*2 "'^Sti 

CllCi2Ci3 



^iie'2 '5*2 «3 ^e3i\ ' 



Mi23 "^I'^b "«2«3 "<^3^i (B16) 

5Cov = except if ^1 = £3 = 0(B17) 

SCov = ^^4^S,^,, 5,,,, 5,3,, (B19) 

Here we do not consider bispectrum coefficients with one 
multipole equal to zero (which amounts to considering the 
power spectrum times the monopole). So the bispectrum 
covariance is diagonal and we find : 



Var(6i23) = ^^ig^^^ X 
JV123 
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